package tech.waterism.algorithm;

import lombok.extern.slf4j.Slf4j;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.interpolation.SplineInterpolator;
import org.apache.commons.math3.analysis.interpolation.UnivariateInterpolator;
import org.springframework.beans.BeanUtils;
import org.springframework.beans.factory.annotation.Autowired;
import org.springframework.stereotype.Component;
import tech.waterism.dto.*;
import tech.waterism.entity.FMParameter;
import tech.waterism.utils.ModelPlatformUtils;

import java.util.ArrayList;
import java.util.Date;
import java.util.stream.Collectors;

@Slf4j
@Component
public class CreatUnitLineMethod {

    private final int unit = 10;
    @Autowired
    private ModelPlatformUtils modelPlatformUtils;

    /**
     * 判断是否水量平衡, 采用水文情报预报规范评估
     *
     * @param runOff    净雨过程 mm
     * @param realFlood 实际洪水过程 m3/s
     * @param area      流域面积 km2
     * @param clen      时间步长 h
     */
    public Boolean checkWaterBalance(double[] runOff, double[] realFlood, double area, double clen) {
        //根据水文情报规范, 径流深的误差的允许范围为实测径流深的20%, 大于20mm, 取20mm, 小于3mm, 取3mm
        double runOffSum = getSum(runOff);
        double criteria = runOffSum * 0.2;
        if (criteria > 20) {
            criteria = 20;
        }
        if (criteria < 3) {
            criteria = 3;
        }
        double runOfferror = getRunOffError(runOff, realFlood, area, clen);
        if (runOfferror > criteria) {
            return new Boolean(false);
        } else {
            return new Boolean(true);
        }
    }

    /**
     * 计算产流误差
     *
     * @param runOff    净雨过程 mm
     * @param realFlood 实际洪水过程 m3/s
     * @param area      流域面积 km2
     * @param clen      时间步长 h
     */
    public double getRunOffError(double[] runOff, double[] realFlood, double area, double clen) {
        //根据水文情报规范, 径流深的误差的允许范围为实测径流深的20%, 大于20mm, 取20mm, 小于3mm, 取3mm
        double runOffSum = getSum(runOff);
        double runOffReal = getSum(realFlood) * clen * 3.6 / area;
        double runOfferror = Math.abs(runOffReal - runOffSum);
        return runOfferror;
    }

    /**
     * 分析法--单位线
     *
     * @param runOff    净雨过程, 第一个数据要求不为0
     * @param realFlood 实测径流过程
     */
    public double[] analysisMethod(double[] runOff, double[] realFlood) throws Exception {
        //        FPInputDto fpInputDto = dto.getFpInputs().get(0);
        //        FPIndexDto fpIndex = dto.getFpIndex();
        //        int size = fpInputDto.getData().size();
        //        double[] runOff = new double[size];
        //        double[] realFlood = new double[size];
        //        for (int i = 0; i < size; i++) {
        //            runOff[i] = fpInputDto.getData().get(i).getRunoffSim();
        //            realFlood[i] = fpIndex.getData().get(i).getQReal();
        //        }
        if (runOff.length > 1 && runOff[1] > runOff[0]) {
            runOff = reverse(runOff);
            realFlood = reverse(realFlood);
        }
        //计算单位线
        double[] unitLine = new double[realFlood.length - 1];
        double unit_h1 = unit / runOff[0];
        for (int i = 0; i < unitLine.length; i++) {
            double tempFlow = 0;
            if (i > 0) {
                int rainLength = i;
                if (rainLength > runOff.length) {
                    rainLength = runOff.length;
                }
                //循环净雨, 累加前面第二个时段净雨开始, 产生的径流
                for (int j = 1; j < rainLength; j++) {
                    tempFlow = tempFlow + unitLine[i - j] * runOff[j] / unit;
                }
            }
            unitLine[i] = (realFlood[i] - tempFlow) * unit_h1;
            if (unitLine[i] < 0) {
                unitLine[i] = 0;
            }
        }
        //倒叙法计算完后，单位线再倒序一下
        if (runOff.length > 1 && runOff[1] > runOff[0]) {
            unitLine = reverse(runOff);
        }
        return formatUnitLine(unitLine);
    }

    /**
     * 试错法做单位线
     *
     * @param fpInputDto 预报方案输入
     * @param clen       计算步长 h
     */
    public double[] tyrErrorMethod(FPInputDto fpInputDto, Integer clen, Date symdh, Date fymdh) throws Exception {
        FMIndexDto fmIndexDto = fpInputDto.getMdcdList().get(1);
        ModelTestDto modelTestDto = new ModelTestDto();
        BeanUtils.copyProperties(fmIndexDto, modelTestDto);
        int size = fpInputDto.getData().size();
        double[] runOff = new double[size];
        double[] realFlood = new double[size];
        for (int i = 0; i < fpInputDto.getData().size(); i++) {
            runOff[i] = fpInputDto.getData().get(i).getRunoffSim();
            realFlood[i] = fpInputDto.getData().get(i).getQReal();
        }
        //计算单位线
        int maxIndex = getMaxValueIndex(runOff);
        double maxValue = getMaxValue(runOff);
        fpInputDto.getData().get(maxIndex).setRunoffSim(new Double(0));
        double[] unitLine_1 = new double[realFlood.length - 1];
        //假设一条单位线UH1
        for (int i = 0; i < unitLine_1.length; i++) {
            unitLine_1[i] = realFlood[i] * unit / maxValue;
        }
        double[] unitLine_2 = new double[realFlood.length - 1];
        double diff = 10;
        int index = 0;
        int totalLineNum = Integer.parseInt(
                modelTestDto.getParameters().stream().filter(i -> i.getNorow().equals(4)).collect(Collectors.toList())
                        .get(0).getParameter().split("//")[0].trim());
        modelTestDto.getParameters().stream().filter(i -> i.getNorow().equals(3)).collect(Collectors.toList()).get(0)
                .setParameter(String.valueOf(totalLineNum + 1));
        int unitLineNo = totalLineNum * 3 + 5;
        FMParameter f1 = new FMParameter().setNorow(unitLineNo).setParameter(String.valueOf(totalLineNum + 1));
        modelTestDto.getParameters().add(f1);
        FMParameter f2 = new FMParameter().setNorow(unitLineNo + 1);
        modelTestDto.getParameters().add(f2);
        FMParameter f3 = new FMParameter().setNorow(unitLineNo + 2);
        modelTestDto.getParameters().add(f3);
        while (diff > 1e-10) {
            //计算去掉后的产流
            //  下面的sim_1为单位线的预报结果, unitLine_1是单位线, unit为单位线的单位, 一般是10mm, ObjectiveFunction.flow_Unit是我自己写的单位线算法
            //  double[] sim_1 = ObjectiveFunction.flow_Unit(runOff_notMax, unitLine_1, unit);
            for (FMParameter parameter : modelTestDto.getParameters()) {
                if (parameter.getNorow() == (unitLineNo + 2)) {
                    StringBuilder sb = new StringBuilder();
                    for (double v : unitLine_1) {
                        sb.append(v);
                        sb.append(" ");
                    }
                    parameter.setParameter(sb.toString());
                }
                if (parameter.getNorow() == (unitLineNo + 1)) {
                    parameter.setParameter(String.valueOf(unitLine_1.length));
                }
            }

            modelPlatformUtils.callModelForecast(modelTestDto, fpInputDto.getData(), clen, null, symdh, fymdh,
                    fpInputDto.getArcd(), fpInputDto.getArtp(), false);
            double[] sim_1 = new double[size];
            for (int i = 0; i < fpInputDto.getData().size(); i++) {
                sim_1[i] = fpInputDto.getData().get(i).getQSim();
            }
            //计算成果与实测的出流量过程相减, 然后计算单位线UH2
            for (int j = 0; j < unitLine_2.length; j++) {
                unitLine_2[j] = (realFlood[j] - sim_1[j]) * unit / maxValue;
                if (unitLine_2[j] < 0) {
                    unitLine_2[j] = 0;
                }
            }
            diff = getDiff(unitLine_1, unitLine_2);
            for (int j = 0; j < unitLine_1.length; j++) {
                unitLine_1[j] = (unitLine_1[j] + unitLine_2[j]) / 2.0;
            }
            //            System.out.println(index++);
            //防止死循环
            if (index > 10000) {
                break;
            }
        }
        return formatUnitLine(unitLine_1);
    }

    /**
     * 单位线时段转换
     *
     * @param sOrgTimeStepTemp 根据单位线计算的S曲线成果
     * @param orgUnitLine      原单位线
     * @param orgTimeStep      原单位线时段长度, h
     * @param newTimeStep      转换后时段长度, h
     */
    public double[] covertUnitLineTimeStep(double[] sOrgTimeStepTemp, double[] orgUnitLine, int orgTimeStep,
                                           int newTimeStep) {
        double[] newUnitLine = null;
        double[] runOff = new double[orgUnitLine.length];
        //默认净雨长度为单位线的长度
        for (int i = 0; i < runOff.length; i++) {
            runOff[i] = 10;
        }
        //S曲线截取
        double[] sOrgTimeStep = new double[orgUnitLine.length + 1];
        double[] simQ_Index = new double[orgUnitLine.length + 1];
        for (int i = 0; i < simQ_Index.length; i++) {
            simQ_Index[i] = i * orgTimeStep;
            sOrgTimeStep[i] = sOrgTimeStepTemp[i];
        }
        int oneLength = (int) simQ_Index[simQ_Index.length - 1] + 1;
        //将S曲线插值成1小时的
        double[] sOneHour = new double[oneLength];
        UnivariateInterpolator insertEntity = new SplineInterpolator();   //单变量内插，样条插值
        UnivariateFunction insertFun = insertEntity.interpolate(simQ_Index, sOrgTimeStep);  //定义单变量函数，插值
        for (int j = 0; j < oneLength; j++) {
            sOneHour[j] = insertFun.value(j);    //插值结果
        }
        //        double[] sNewTimeStep = new double[oneLength];
        ArrayList<Double> temp = new ArrayList<Double>();
        temp.add(new Double(0));
        for (int j = 0; j < oneLength; j++) {
            double sNewTimeStep = 0;
            if (j >= newTimeStep) {
                //s(t) - s(t - newTimeStep)
                sNewTimeStep = sOneHour[j] - sOneHour[j - newTimeStep];
                //q(t) = orgT/newT*(s(t) - s(t - newTimeStep))
                sNewTimeStep = orgTimeStep * sNewTimeStep / newTimeStep;
            }
            //获取对应时段并且大于0的单位线数据
            if (j % newTimeStep == 0) {
                if (sNewTimeStep > 0) {
                    temp.add(new Double(sNewTimeStep));
                }
            }
        }
        temp.add(new Double(0));
        //最后将List转换为数据，即为单位线
        newUnitLine = new double[temp.size()];
        for (int j = 0; j < newUnitLine.length; j++) {
            newUnitLine[j] = temp.get(j).doubleValue();
        }
        return newUnitLine;
    }

    /**
     * 单位线水量平衡调整
     *
     * @param unitLine 单位线数据
     * @param area     流域面积，平方千米
     * @param timeStep 时间步长，小时
     */
    public double[] unitLineWaterBalance(double[] unitLine, double area, int timeStep) {
        //单位线的Q
        double sumQ = getSum(unitLine);
        //根据单位线的径流深计算的Q
        double standardQ = unit * area / (3.6 * timeStep);
        //计算单位线10mm对应的Q与单位线计算的Q的比例d
        double diff = standardQ / sumQ;
        //等比例修正单位线
        for (int i = 0; i < unitLine.length; i++) {
            unitLine[i] = unitLine[i] * diff;
        }
        return unitLine;
    }

    /**
     * 计算单位线的水量
     *
     * @param unitLine 单位线
     * @param timeStep 时间步长 h
     * @return 单位线下包的水量，万方
     */
    public double getUnitLineVolume(double[] unitLine, int timeStep) {
        //单位线的Q
        double sumQ = getSum(unitLine);
        double volume = 3.6 * sumQ * timeStep;
        return volume;
    }

    private double getDiff(double[] x, double[] y) {
        double diff = 0;
        for (int i = 0; i < x.length; i++) {
            diff = diff + Math.abs(x[i] - y[i]);
        }
        return diff;
    }

    private double getMaxValue(double[] q) { // 取最大值（float）
        if (q == null)
            return 0;
        double max = 0;
        for (int i = 0; i < q.length; i++) {
            if (q[i] > max)
                max = q[i];
        }
        return max;
    }

    private int getMaxValueIndex(double[] q) {
        if (q == null)
            return 0;
        double max = 0;
        int maxIndex = 0;
        for (int i = 0; i < q.length; i++) {
            if (q[i] > max) {
                max = q[i];
                maxIndex = i;
            }
        }
        return maxIndex;
    }

    private double getSum(double[] q) {
        if (q == null)
            return 0;
        double sum = 0;
        for (int i = 0; i < q.length; i++) {
            sum = sum + q[i];
        }
        return sum;
    }

    private double[] reverse(double[] array) {
        double[] newArray = new double[array.length];
        int index = 0;
        for (int i = array.length - 1; i >= 0; i--) {
            newArray[index] = array[i];
            index++;
        }
        return newArray;
    }

    private double[] formatUnitLine(double[] array) {
        ArrayList<Double> list = new ArrayList<>();
        list.add(new Double(0));
        for (int i = 0; i < array.length; i++) {
            if (array[i] >= 1) {
                list.add(array[i]);
            }
        }
        list.add(new Double(0));
        double[] result = new double[list.size()];
        for (int i = 0; i < result.length; i++) {
            result[i] = list.get(i);
        }
        return result;
    }

    /**
     * Nash and Sutcliffe, 确定性系数
     *
     * @param real 实际值
     * @param sim  模拟值
     */
    private double NSFun(double[] real, double[] sim) {
        double NSvalue = 0;
        double real_avg = getSum(real) / (double) real.length;
        double sum1 = 0;
        double sum2 = 0;
        for (int i = 0; i < sim.length; i++) {
            sum1 += Math.pow(real[i] - sim[i], 2);
            sum2 += Math.pow(real[i] - real_avg, 2);
        }
        NSvalue = 1 - sum1 / sum2;
        return NSvalue;
    }
}
